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Abstract The Bayesian approach to inverse problems is of paramount impor- 
tance in quantifying uncertainty about the input to and the state of a system 
of interest given noisy observations. Herein we consider the forward problem 
of the forced 2D Navier Stokes equation. The inverse problem is inference of 
the forcing, and possibly the initial condition, given noisy observations of the 
velocity field. We place a prior on the forcing which is in the form of a spa- 
tially correlated temporally white Gaussian process, and formulate the inverse 
problc^ni for the posterior distribution. Given appropriate spatial regTilarity 
conditions, we show that the solution is a continuous function of the forcing. 
Hence, for appropriately chosen spatial regularity in the prior, the posterior 
distribution on the forcing is absolutely continuous with respect to the prior 
and is hence well-defined. Furthermore, the posterior distribution is a continu- 
ous function of the data. We complement this theoretical result with numerical 
simulation of the posterior distribution. 



1 Introduction 

The Bayesian approach to inverse problems has grown in popularity signif- 
icantly over the last decade, driven by algorithmic innovation and steadily 
increasing computer power [7]. Recently there have been systematic develop- 
nu^uts of tlu- tlKX)ry of Bayesian inversion on function space [8,9,2,15,10,11] 
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and this has led to new sampUng algorithms which perform well under mesh- 
refinement [1, 17]. In this paper we add to this growing interest in the Bayesian 
formulation of inversion, in the context of a specific PDE inverse problem, mo- 
tivated by geophysical applications such as atmospheric/occanographic data 
assimilation, and demonstrate that fully Bayesian probing of the posterior 
distribution is feasible. 

The primary goal of the paper is to demonstrate that the Bayesian for- 
mulation of inversion for the forced Navier-Stokes equation, introduced in [2] , 
can be extended to the case of white noise forcing; the paper [2] assumed 
an Ornstcin-Uhlenbeck structure in time for the forcing, and hence did not 
include the white noise case. It is technically demanding to extend to the 
case of white noise forcing, but it is also of practical interest. This practical 
importance stems from the fact that the Bayesian formulation of problems 
with white noise forcing corresponds to a statistical version of the continu- 
ous time weak constraint 4DVAR methodology [18]. The 4DVAR approach to 
data assimilation currently gives the most accurate global short term weather 
forecasts available [13] and this is arguably the case because, unlike ensemble 
filters which form the major competitor, 4DVAR has a rigorous statistical in- 
terpretation as a MAP estimator. It is therefore of interest to seek to embed 
our understanding of such methods in a broader Bayesian context. 

The key tools required in applying the function space Bayesian approach 
in [15] are the proof of continuity of the forward map from the function space 
of the unknowns to the data space, together with estimates of the dependence 
of the forward map upon its point of application, sufficient to show certain in- 
tegrability properties with respect to the prior. This program is carried out for 
the 2D Navier-Stokes equation with Ornstein-Uhlenbeck priors on the forcing 
in the paper [2]. However to use priors which are white in time adds fur- 
ther complications since it is necessary to study the stochastically forced 2D 
Navier-Stokes equation and to establish continuity of the solution with respect 
to small changes in the Brownian motion which defines the stochastic forcing. 
We do this by employing the solution concept introduced by Flandoli in [4], 
and using probabilistic estimates on the solution derived by Mattingly in [14] . 
In section 2 we describe the relevant theory of the forward problem, employing 
the setting of Flandoli. In section 3 we build on this theory, using the estimates 
of Mattingly to verify the conditions in [15], resulting in a well-posed Bayesian 
inverse problem for which the posterior is Lipschitz in the data with respect 
to Hellinger metric. Section 4 extends this to include making inference about 
the initial condition as well as the forcing. Finally, in section 5, we present 
numerical results which demonstrate feasibility of sampling from the posterior 
on white noise forces. 



2 Forward Problem 

In this section we study the forward problem of the Navier-Stokes equation 
driven by white noise. Subsection 2.1 describes the forward problem, the 
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Navier-Stokcs equation, and rewrites it as an ordinary differential equation 
in a Hilbert space. In subsection 2.2 we define the functional setting used 
throughout the paper. Subsection 2.3 highlights the solution concept that we 
use, leading in subsection 2.4 to proof of the key fact that the solution of 
the Navier-Stokes equation is continuous as a function of the rough driving of 
interest. 



2.1 Overview 

Let e be a bounded domain with smooth boundary. We consider in D 
the Navier-Stokes equation 

dtU-vAu + u-Vu = f -Vp, (x, t) e D X (0, oo) (1) 
V-u = 0, (a;,<) e D X (0,00) 
M = 0, (a;,t) e 9£) X (0,00), 
u = uq, {x,t) € D X {0}. 

We assume that the initial condition uq and the forcing f{-,t) are divergence- 
free. We denote by V the space of all divergcncc-frce smooth functions from 
D to with compact support, by H the closure of V in and by 

the closure of V in {H\D)f. Let = (H^ (D))'^ f]m\ The initial condition 
Mo is assumed to be in H. We define the linear Stokes' operator A : ^ H 
by Au = —Au noting that the assumption of compact support means that 
Dirichlet boundary condition are imposed on the Stokes' operator Since A is 
selfadjoint, A possesses eigenvalues < Ai < A2 < . . . with the corresponding 
eigenvectors ei, 62, . . . G H^. 

We denote by (•, •) the inner product in H, and extend to a dual pairing 
on X We then define the bilinear form S : x ^ H'^ 

{B{u,v),z)— / z{x) ■ {u{x) ■ V)v{x)dx 
Jd 

which must hold for all z ^V}. From the incompressibility condition we have, 
for all z G Hi, 

{B{u,v),z) =-{B{u,z),v). (2) 

By projecting problem (1) into H we may write it as an ordinary differential 
equation in the form 

du{t) = -vAudt - B{u, u)dt + dW{t), u(0) = uo G H, (3) 

where dW{t) is the projection of the forcing /(x, t)dt into H. We will define the 
solution of this equation pathwise, for suitable W, not necessarily differentiable 
in time. 



If, instead, we consider the problem on the torus T^, all the results below hold. 
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2.2 Function Spaces 

For any s > we define c H to be the Hilbert space of functions u = 
Yl'kLi Uh^k € H such that 

oo 

Yl ^k^l < oo; 

k=l 

we note that the for j e {0, 1, 2} coincide with the preceding definitions of 
these spaces. The space is endowed with the inner product 

oo 

fc=l 

for u = X]fe°=i '^k^k, V = X]fc°=i ^feGfe G EI- Wc denote by V the particular choice 
s = \ + t, namely El5+'^, for e > 0. In what follows we will be particularly 
interested in continuity of the mapping from the forcing W into linear func- 
tional of the solution of (3). To this end it is helpful to define the Banach 
space X := C([0,T];V) with the norm 

\\W\U= sup \\W{t)h- 
te{o,T) 



2.3 Solution Concept 

In what follows we define a solution concept for equation (3) in the case where 
W is continuous, but not necessarily differentiable, in time. We always assume 
that W{0) = 0. The approach is based on that in Flandoli [4]. In this approach 
we first make sense of the equation 

dz{t) = -uAzdt + dW{t), 2(0) = 0eH. (4) 

Then we write u = z + v and notice that v solves the equation 

dv{t) = -vAvdt- B{z + v,z + v)dt, d(0) = e H. (5) 

We make sense of this equation for v and then deduce that u = z + v is a. weak 
solution (defined below) of (3). When we wish to emphasize the dependence 
of w on (and similarly for z and v) we write u{t; W). 

For each e X, we define the weak solution u{-, W) e C([0, T]; H) f] ^^([o, T]; H^/^) 
of (3) as a function that satisfies 

{u{t),cj)}+u [ {u{s),Acl))ds- I {B{u{s),4>),u{s))dx = {uQ,4>) + {W{t),4>), 
Jo Jo 

(6) 

for all ^ gIP and all t G (0, T); note the integration by parts on the Stokes' 
operator and the use of (2) to derive this identity from (3). 
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Application of the variation of constants formula to (4) , together with an 
integration by parts, yields 



z{t) 



-uA(t—s 



= W{t) - [ uAe-'''^'^*-'^W{s)ds. 
Jo 

We now justify this formula for z. We define Wk = {W, e^), that is 

oo 

W{t) :=^w;fc(t)efe ex. 

k=l 

We then define 

oo 

z{t;W) = W{t) Wk{s)iyXke^'-'^^-''^>'Usyk, 



(7) 



(8) 



fe=i 



noting that this agrees with the expression for z{t) above. We then have the 
following: 

Lemma 1 For each W eX, the function z = W) E C([0, T]; Hi/^). 

Proof We first show that for each t, z{t;W) as defined in (8) belongs to 
H^/^. Fixing an integer M > 0, using inequality a^~^/^e~" < c for all a > 
for an appropriate constant c, we have 



fe=i V-^o 
Therefore, 

V / uXkt 



Wk 



,(f-s)(-i^Afc) 



Wk{s)ekds 



Hi/2 



< 



< 



■e/2 

M 



Xk'^\'Wk{s)\ekds 



^A^/Vfc(s)|efe 



Hi/2 

ds 



< niax ||H^(s)||ui/2+e / 
which is uniformly bounded for all M. Therefore, 



(i-s)W2 



Hi/2 



ds, 



It follows from (8) that, since € X, for each t, z{t; W) e H^/^ as required. 
Furthermore, for all t e (0, T) 



|^(t;W^)||Hi/2 <c||W^||x. 



(9) 
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Now we turn to the continuity in time. Arguing similarly, we have that 



/o (i-s)W2 

O 

Wk{s)ek 



k=M 



Wk{s)ek 

k=M 
1 \ 1/9 

ds 



ds 



Hl/2 + e 



Hi/2+« 



for all p,q> Q such that 1/^+1/^=1. Prom the Lebesgue dominated con- 
vergence theorem, 

9 



lim 



Jo 



Wk{s)ek 



k=M 



ds = 0; 



Hl/2+€ 



and when p sufficiently close to 1, 



-e/2)p 



ds 



is finite. We then deduce that 



lim 

M—^oo 



k=M \"^° 



1/2 



uniformly for all 

Fixing t G (0, T) we show that 

\im\\z{t■,W)-z{t'■,W)\\^^,.=0. 

Wc have 

W) - z(t'- W)Uu^ < \\W{t) - VF(t')llHV2 + 

^ ^jl'"k{s)y\ke^'-'^^-''^''^ds- j'^ Wk{s)u\ke^'-'^^-"^''^d^ek 



+ 



Hi/2 



£ ( f Wk{s)v\ke('-'^^-''^''^ds\ek 
k=M\-^o J 



J2 I I Wk{s)uXke^'-'^^-'^''Us]ek 



k=M 



Hi/2 



For (5 > 0, when M is sufficiently large, the argument above shows that 

£ ( f Wk{s)v\ke('-'^(-''^''Us\ek + 
Furthermore, when — i| is sufficiently small, 

E 



k=M\''0 J 



<6/3. 



M-1 / ,t ,t' \ 



fe=i 



Jo 



< 6/3. 
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Finally, since W gX, for \t' — t\ is sufficiently small we have 

\\W{t)-W{t')\\j,^,2<d/3 

Thus when \t' - t\ is sufficiently small, - z{t';W)\\^i/2 < 6. The 

conclusion follows. □ 
Having established regularity, we now show that z is indeed a weak solution 
of (4). 

Lemma 2 For each <p gW^, z{t) = z{t; W) satisfies 

{z{t), cj>)+u f\z{s),Acl))ds = {W{t), 4,). 
Jo 

Proof It is sufficient to show this for ^ = efe. We have 

/ {zis),Aek)ds= f {W{s),Aek)ds- I f Wk{T)vXle^'-''^<^-''^''UTds 
Jo Jo Jo Jo 

= Afe^ Wk{s)ds - uXl Wk{T)(^j e^'-'^^'^-^^'Us^dT 

= Xkf Wk{s)ds-Xkf Wk{T)dT + Xk f M;fe(T)e(*-^)(-"^'=)dT 
Jo Jo Jo 

Jo 

On the other hand, 

{z{t),ek) = {W{t),ek) - vXk f Wk{s)e('-'^(-''^>'Us. 

Jo 

The result then follows. □ 

Wc now turn to the following result, which concerns v and is established 
on page 416 of [4], given the properties of z{-;W) estabhshed in the preceding 
two lemmas. 

Lemma 3 For each e X, problem (5) has a unique solution v G C(0, T; H) f] L'^{0, T 

We then have the following existence and uniqueness result for the Navier- 
Stokes equation (3), more precisely for the weak form (6), driven by rough 
additive forcing [4]: 

Proposition 1 For each W G X, problem (6) has a unique solution u G 
C(0,T;H)ni^(0,T;Hi/2) such that u - z G L^{0,T;M^). 

Proof A solution u for (6) can be taken as 

u{t;W) ^ z{t:W) +v{t:W). (10) 

Assume that u{t;W) is another solution of (6). Then v{t;W) = u{t;W) — 
z{t;W) is a solution in C(0,T;H) 0^^(0,7; H^) of (5). However, (5) has a 
unique solution in C(0, T; H) fl L'^{0, T; H^). Thus v = v. □ 
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2.4 Continuity of the Forward Map 



The purpose of this subsection is to establish continuity of the forward map 
from W into the weak solution u of (3), as defined in (6), at time t > 0. 

Theorem 1 For each t > 0, the solution u{t; •) of (3) is a continuous map 
from X into H. 

Proof Wc consider equation (3) with driving W G X given by (7) and by 
W gX defined by 

OO 

W\s) = J2^k'{s)ekeX. 
fe=i 

We will prove that, for W, W from a bounded set in X, there is c = C{T) > 0, 
such that 

sup \\z{t; W) - z{t; W')U^J. < c\\W - W'\\x (11) 
te(o,T) 

and, for each t G (0,T), 

\\v{t;W)-v{t;W')\g<c sup Us;W) - z{s;W')\\l.^j,y (12) 

sG(0,T) 

This suflBces to prove the desired result since Sobolev embedding yields, from 
(12), 

\\v{t;W)-v{t;W')\\l<c sup \\z{s;W) - zis;W')\\\ . (13) 

se(o,T) 



'H2 



Since u = z + v we deduce from (11) and (13) that u as a map from X to . 
is continuous. 

To prove (11) we note that 



\\z{t;W)-z{t;W')\\^. < \\W{t)-W' m^.+ 



so that 



sup \\z{t; W) - z{t; W')\\^, < \\W - W 



vAe-'^^'-''^ {W{s)-W'{s))ds 



H2 



te(o,T) 



- sup 
te(o,T 



f i/^e-"^(*-^) {W{s) - W'{s))ds 
Jo 

Thus it sufiices to consider the last term on the right hand side. We have 

f Ae-''^^'''\w{s)-W'{s))ds 
Jo 



H2 



^1 ^0 



(t-s){-v\k)(.,,i 



{w'kis) - Wk{s))ekds 



Hi/2 
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-S)(-L/Afc)(- / 



(w)fc(s) - Wk{s))ds 



< 



< 



where we have used the fact that ^/^e " < c for all a > for an appropriate 
constant c. From this, we deduce that 



^g-.A(t-.) (w{s)-W'{s))ds 

OO 



H2 



< 



< 



fe=l 





< c sup ||VF'(s) - W{s) 
se{o,T) 



Hi/2 



Hi/2 



< 



ds sup 



fc=i 



Therefore (11) holds. 

Wc now prove (12). Wc will use the following estimate for the solution 
V of (5) which is proved in Flandoli [4], page 412, by means of a Gronwall 
argument: 

sup \\v{s)\\m+ r \\v{s)rs.<C{T, sup \\z{s)U.^D)). (14) 
se(o,T) Jo se(o,T) 



Wc show that the map C([0,r];i''(D)) 3 z{-,W) ^ v{-,W) G H is continuous. 
For W and W in X, define v = v{t;W),v' = v{t;W'),z = z{t;W),z' = 
z{t; W), e = V — v' and 6 = z — z'. Then we have 



de 
di 



+ vAe -\-B{y-\-z,v-\-z)-B{y' + z' , v' + z') = 0. 



(15) 



From this, we have 



IdWef 



2 dt 
-{B[v + z,v + z),e) 

+{B{v' + z',v' + z'),e). 
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Prom (2) we obtain 



lrf||e||^ 



■i^l|e||Hi = 



2 dt 
+{B{v + z,e),v + z) 

-{B{v' + z',c^,v' + z') 

= {B{v + z,e),v + z - v' - z') 

-{B{v' + z' -V- z, e),v' + z') 

= {B{v + z,e),e + 6) 
+{B{e + S, e),v' + z') 
< i\\4LHD) + \\S\\lhd)){\\v\\lhd) + \\z\\L-iD) + \\v'\\l^(d) + \\z'\\L^(D))\\e\\m^- 
Using the interpolation incquahty 

\\e\\L^iD)<c\\ed\e\\'J', 

we get 

2 dt +^"11^1111^^ 

c(l|e|||(^||e||J/^ + ||<5|U4(r,)||e||Hi) ■ (^Ml^d) + \\v'\\l^(d) + Ml^d) + W- 
Prom the Young inequality, we have 

I ^ + H|e||^. < u\\e\g, + c\\e\g ■ h + c||<5||i.(^) • h- (16) 
where we have defined 



From Gronwall's inequality, we have 

< cf^ {ers^^(^')'i^y5{s)\\l.^^)h{s)ds. (17) 
Prom the interpolation inequality 

\\v{s'; W)|U.(Z5) < c\Hs'- W)\\]ihv{s'; W)\\U\ 

we have that 

f \\v{s';W)\\i.^j,^ds' <c sup \\v{s';W)\\l f \\v{s';W)\\lM, 
Jo s'e(o,T) Jo 

which is bounded uniformly when W belongs to a bounded subset of X due to 
(14) and (9). Therefore 



||e(t)||^<c sup ||5(s)||i4(B). 

0<s<T 



□ 
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3 Bayesian Inverse Problems With Model Error 

In this section we formulate the inverse problem of determining the forcing to 

equation (3) from knowledge of the velocity field; more specifically we formu- 
late the Bayesian inverse problem of determining the driving Brownian motion 
W from noisy pointwise observations of the velocity field. 

We set-up the likelihood in subsection 3.1. Then, in subsection 3.2, we 
describe the prior on the forcing which is a Gaussian white-in-time process with 
spatial correlations, and hence a spatially correlated Brownian motion prior 
on W. This leads, in subsection 3.3, to a well-defined posterior distribution, 
absolutely continuous with respect to the prior, and Lipschitz in the Hellinger 
metric with respect to the data. 

3.1 Likelihood 

Fix a collection of times tj G {0,T), j = 1, - ■ ■ , J. Let £ be a collection of K 
linear functionals on H. We assume that we observe, for each j, £{u{-, tj; W)) 
plus a draw from a centered Gaussian noise ^? with variance a, i.e. 



where (5, <E M"^^ and : X ^ R"^^ . We assume that the observational noise 
■d is a, draw from the Gaussian A^(0, U) on R''^ . 

In the following we will define a prior measure ponW and then determine 
the conditional probability measure = P{W\S). We will then show that 

is absolutely continuous with respect to p and that the Radon-Nikodym 
derivative between the measures is given by 



5j=i{u{;tj;W))+^j, 



(18) 



is known to us. Concatenating the data we obtain 



(19) 




(20) 



where 




(21) 



The right hand side of (20) is the likelihood of the data S. 



3.2 Prior 



We construct our prior on the time-integral of the forcing, namely W. Let Q 
be a linear operator from the Hilbert space into itself with eigenvectors 
Cfe and eigenvalues cr^ for A; = 1, 2, We make the following assumption 



12 



Viet Ha Hoang et al. 



Assumption 1 There is an e> such that the coefficients {ak} satisfy 

oo 
fe=l 

As 

OO OO _^ 

fe=l /c=l 

Q is a trace class operator in 

We assume that our prior is the Q- Wiener process W with values in H3+^ 
where W{si) — W{s2) is Gaussian in IHl2+^ with covariance (si — S2)Q and 
mean 0. In the mean square norm, this process can be written as 

oo 

w{t)^j2'^>'^km), (22) 
fe=i 

where /3fe(t) are pair- wise independent Brownian motions (sec Da Prato and 
Zabczyk [3], Chapter 4). We define by p the measure generated by this Q- 
Wiener process on X. 

Remark 1 We now show that under Assumption 1, p almost surely, solution 
u of (6) defined in (10) equals the unique progressively measurable solution in 
C([0,T];H)n-£'^([0,T];Hi) constructed in Flandofi [4] when the noise W is 
sufficiently spatially regular. This enables to employ the boundedness of the 
second moment of the energy E^[||u(-, T4^)||^], established in Mattingly [14], 
which we need later. 

For the infinite dimensional Brownian motion W defined in (22) where 

oo 

fc=l 

for some /3o > 0, Flandoli [4] employs the Ornstein-Uhlenbeck process 

Zoc{t)= f e-("^+")(*-«)rfW^(s) (23) 

J —oo 

where a is a constant in order to define a solution of (6). Note that if /3o > 5 
then Assumption 1 is satisfied. With respect to the probability space (i7, J^t, P), 
the expectation E|jzQ,(t)||giy2_|_2^ is finite for /? < Pq. Thus almost surely with 
respect to (i7, J"t,P), Za{f) is sufficiently regular so that problem (5) with the 
initial condition v{Q; W) — uq — ^q(O) is well posed. The stochastic solution to 
the problem (1) is defined as 

u{;t;W) = Zo,{t;W)+v{t;W). (24) 

When^o > \,^\zo,{t)\\'^, is finite so u{-,t;W) & C{[Q,T]-W)[^L'^{[0,T]-m^). 
Flandoli [4] leaves open the question of the uniqueness of a generalized solution 
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to (6) in C([0,T]:H) ni^([0,T];HV2). However, there is a unique solution in 
C([0,r];H)ni'([0,r];Hi). 

Almost surely with respect to the probability measure p, solution u of (6) 
constructed in (10) equals the solution constructed by Flandoli [4] in (24). To 
see this, note that the stochastic integral 

t 

e-''^(*-")dW(s) (25) 
can be written in the integration by parts form (8). Therefore, with respect to 



W[\\z{t)\\ 




which is finite. Therefore p almost surely, z{t) e L^(0, T;!!^). Thus uif/, W) G 
C(0,T;H)ni^(0.r:Hi). We can then argue that p almost siuely, the solu- 
tion u constructed in (10) equals Flandoli's solution in (24) which we de- 
note by Ua (even though it does not depend on a) as follows. As it^ G 
C([0, T];H)f] i2([0, T]; Hi), v^{t; W) = u„(i; W)-zit; W) G C([0, T];H)f] L^{[0, T]; H^) 
and satisfies (5). As for each W, (5) has a unique solution in C([0, T];H)f] L^{[0, T]; H^), 
so Va{t; W) = v{t; W). Thus almost surely, the Flandoli [4] solution equals the 
solution u in (10). This is also the argument to show that (3) has a unique 
solution in C( [0, T] ; H) fl ( [0, T] ; Hi ) . 



3.3 Posterior 

Theorem 2 The conditional measure F{W\5) = p^ is absolutely continuous 
with respect to the prior measure p with the Radon-Nikodym derivative being 
given by (20) . Furthermore, there is a constant c so that 

dHen{p\p'')<c\S-5'\. 

Proof From Corollary 2.1 of Cotter et al. [2], it suffices to show that the 
mapping (J : X — > M'^^^ is continuous, in order to deduce that p"* <§; p and that 
the Randon-Nikodym derivative (20) holds. This follows from Theorem 1. The 
Lipschitz continuity of the posterior in the Hellingcr metric follows similarly to 
the proof of Theorem 2.5 in [2] and we only sketch the proof. As in that proof 
we observe that it is sufficient to consider observations at one time to > 0. We 
define 

Z{6) := [ exp{-^{W;S))dp{W). 
Jx 
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Mattingly [14] shows that for each t, E''(||w(-, t; W)||g) is bounded. Fixing 
a large constant M, the probabihty that W^)||h < M is larger than 

1 - c/M > 1/2. For such a path W, 

$iW:S)<ci\S\ + \\e\\M-iM). 

From this, we deduce that Z{6) > 0. Next, we have that 

\Z{5)- Z{d')\ < [ \${W;5)-${W]5')\dp{W) 



< c 



/ {\5\ + |<5'| + 2KKio; W))\)\5 - 5'\dp{W) 
Jx 



< c\S-6'\. 

We then have 



2dnen{p',p''f < ^ (^Z{6)-'/^eM-lHW;S)) - ^(5')"'/' exp(-i#(l^'; 5'))) dp{W) 



< h+h, 

where 



and 



= 2\Z{5)-^'^ - Z(y)-V2|2 j exp(-^(W; 5'))dp{W). 
From this, we deduce that dHeii(p'^) P*^ ) < c\5 — 5'\. 



□ 



4 Inferring The Initial Condition 

In the previous section we discussed the problem of inferring the forcing from 
the velocity field. In practical applications it is of interest to infer the ini- 
tial condition, which corresponds to a Baycsian interpretation of 4DVAR, or 
the initial condition and the forcing, which corresponds to a Bayesian inter- 
pretation of weak constraint 4DVAR. Thus we consider the Bayesian inverse 
problem for inferring the initial condition uq and the white noise forcing deter- 
mined by the Brownian driver W . Let be a Gaussian measure on the space 
H and let /x = p be the prior probability measure on the space ?^ = H x X. 
We denote the solution u of (3) by u{x, t; uq, W). 

We outline what is required to extend the analysis of the previous two 
sections to the case of inferring both initial condition and driving Brownian 
motion. We simplify the presentation by assuming observation at only one 
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time to > although this is easily relaxed. Given that at to € (0,T), the noisy 
observation S of £{u{-, to;uo, W) is given by 

S = £{u{;to;uo,W))+^ (26) 

where 'd ~ iV(0, S). Letting 

#(«o, W;S) = ^\S- i{u{; to; «o, W))\l, (27) 
we aim to show that the conditional probability yu* = P('Uoj W\5) is given by 

ocexp{-<P{uo,W;5)). (28) 

We have the following result. 

Theorem 3 The conditional probability measure ji^ = V(uo,W\5) is abso- 
lutely continuous with respect to the prior probability measure ^ with the Radon- 
Nilodym derivative give by (28) . Further, there is a constant c such that 

dHeu{lJ'\/)<c\5-5'\. 

Proof We only sketch the proof because details are similar to those in the 
previous section. The key issue is establishing continuity of the forward map 
with respect to initial condition and driving Brownian motion. We show that 
u{- J/,uo^W) is a continuous map from % to H. For g X and Uo € H, we 
consider the following equation: 

dv 

— +Av + B{v + z,v + z) = 0, v{0) = uo. (29) 

We denote the solution by v{t) = v{t\ uq, W) to emphasize the dependence on 
initial condition and forcing which is important here. For (uq, W) G Ji and 
(uq, W) e H, from (16) and Gronwall's inequality, we deduce that 

\\v{t:uo,W)-v{t:u'„W')\\l < \\uo - uX^eJo ^^(^'))'^^' 
+cf^ (e/^^(^')''^' ■ \\z{s;W)- z{s;W')\\l.(^yh{s)ys. 
We then deduce that 

\\v{t; uo, W) - v{t; u'o, W')\g < c\\uo - u'o\\l + c sup_ \\z{s; W) - z{s; W')\\l.^j,^ 



0<s<T 

.1 l|2 



< c||uo-Uo||h + c sup \\W{t)-W'{t)\\^u2+.. 

te(o,T) 

This gives the desired continuity of the forward map. 

For the Lipschitz dependency of the Hellinger distance of /x'^ on S, we use 
the result of Mattingly [14] which shows that, for each initial condition Uq, 

EPi\\u{t;uo,W)\\l,) < A_ +e-2^^^*(||«o||^- ^), 

where So = J2kLi '^k- Therefore E'*(||u(t; uo, W) ||^) is bounded. The remainder 
of the proof follows as in Theorem 2. □ 
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5 Numerical Results 

The purpose of this section is to demonstrate that the Bayesian formulation 
of the inverse problem described in this paper forms the basis for practical 
numerical inversion. In particular we show that it is possible to recover white 
noise forcing of the Navier-Stokes equation from linear functional of the ve- 
locity field. In subsection 5.1 we describe the numerical method used for the 
forward problem. In subsection 5.2 we describe the inverse problem and the 
Metropolis-Hastings MCMC method used to probe the posterior. Subsection 
5.3 describes the numerical results. 



5.1 Forward Problem: Numerical Discretization 

All our numerical results are computed using a viscosity of = 0.1 and on the 
periodic domain. We work on the time interval t G [0,0.1]. We use M = 32^ 
divergence free Fourier basis functions for a spectral Galerkin spatial approxi- 
mation, and employ a time-step 5t = 0.01 in a Taylor time-approximation [6]. 
The number of basis functions and time-step lead to a fully-resolved numerical 
simulation at this value of v. 



5.2 Inverse Problem: Metropolis Hastings MCMC 

We consider the inverse problem of finding the driving Brownian motion. As 
a prior we take a centered Brownian motion in time with spatial covariance 
TT^A"^; thus the space-time covariance of the process is Co := Tr^A~'^ (g) A^^, 
where At is the Laplacian in time with fixed homogeneous Dirichlet condition 
ait = and homogeneous Neumann condition ait = T. It is straightforward to 
draw samples from this Gaussian measure, using the fact that A is diagonalized 
in the spectral basis. Note that iiW ^ p, then W G C(0, T; W) almost surely 
for all s < 1; in particular W gX. Thus p{X) = 1 as required. The likelihood 
is defined by making observations of the velocity field at every point on the 
32^ grid implied by the spectral method, at every time t = n6t, n = 1, • • • 10. 
The observational noise standard deviation is taken to be 7 = 3.2. 

To sample from the posterior distribution we employ a Metropolis-Hastings 
MCMC method. Furthermore, to ensure mesh-independent convergence prop- 
erties, we use a method which is well-defined in Hilbert space [1]. Metropolis- 
Hastings methods proceed by constructing a Markov kernel V which satisfies 
detailed balance with respect to the measure which we wish to sample: 

p\du)r{u, dv) = p\dv)V{v, du), V u, u e X. (30) 

Integrating with respect to u, one can see that detailed balance implies p^V = 
p*. Metropolis-Hastings methods [5,16] prescribe an accept-reject move based 
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on proposals from another Markov kernel Q, in order to define a kernel V 
which satisfies detailed balance. If we define the measures 



u^{du,dv) — Q{v,du)p°{dv) oc Q{v,du)expl—${v;S)]p{dv). 

then, providc;(l i^-^ <^ v, the Metropolis-Hastings method is defined as follows. 
Given current state u„, a proposal is drawn u* ~ Q(wn, •)> ^"^^ then accepted 
with probability 



Q!(m„,u*) = min |l, -^(u„,u*)| . (32) 

The resulting chain is denoted by P. If the proposal Q preserves the prior, so 
that pQ = p, then a short calculation reveals that 

a{un,u*) = min |l,exp^^(M„; 5) - (l>{u*;S)j | ; (33) 

thus the acceptance probability is determined by the change in the likelihood 
in moving from current to proposed state. We use the following pCN proposal 
[1] which is reversible with respect to the Gaussian prior 

Q(«„, •) = ^i^/l^un, l3^Co) (34) 

and thus results in the acceptance probability (33). Variants on this algorithm, 

which propose differently in different Fourier components, arc described in [12], 
and can make substantial speedups in the Markov chain convergence. However 
for the examples considered here the basic form of the method suffices. 



5.3 Results and Discussion 

The true driving Brownian motion , underlying the data in the likelihood, 
is constructed as a draw from the prior p. We then compute the corresponding 
true trajectory u'^{t) = u{t; W^). We use the pCN scheme (33), (34) to sample 
W from the posterior distribution p^. The true initial and final conditions are 
plotted in Figure 1, top two panels, for the vorticity field w; the bottom two 
panels of Figure 1 show the posterior mean of the same quantities and indicate 
that the data is fairly informative, since they closely resemble the truth. The 
true trajectory, together with the posterior mean and one standard deviation 
interval around the mean, are plotted in Figure 2, for the wavenumbers (0, 1) 
and (4, 4) , and for both the driving Brownian motion W and the velocity field 
u. This figure indicates that the data is very informative about the (0, 1), but 
less so concerning the (4, 4) mode for which the mean and standard deviation 
exhibit behaviour similar to that under the prior. 
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Fig. 1 The truth (top) and expected value (bottom) of u{t;W), with t = 0.01 (left) and 
t = 0.1 (right). 



References 

1. S. Cotter, G. Roberts, A. Stuart, and D. White. MCMC methods for functions; modi- 
fying old algorithms to make them faster. Arxiv preprint arXiv: 1202.0709, 2012. 

2. S. L. Cotter, M. Dashti, J. C. Robinson, and A. M. Stuart. Bayesian inverse problems 
for functions and applications to fluid mechanics. Inverse Problems, 25:115008, 2009. 



Determining White Noise Forcing of the Navier Stocks Equation 



19 





Fig. 2 The trajectories Uk{t] W) (top) and (bottom), with k = (0, 1) (left) and k = (4, 4) 
(right). Shown are expected values and standard deviation intervals as well as true values. 



3. Guiseppe Da Prato and Jerzy Zabczyk. Stochastic Equations in Infinite Dimensions. 
Cambridge University Press, 2008. 

4. F. Flandoli. Dissipative and invariant measures for stochastic Navicr-Stokcs equations. 
NODEA, 1:403-423, 1994. 

5. W.K. Hastings. Monte carlo sampling methods using markov chains and their applica- 
tions. Biometrika, 57:97-109, 1970. 



20 



Viet Ha Hoang et al. 



6. A. Jcntzcn and P. Klocdcn. Taylor expansions of solutions of stochastic partial differ- 
ential equations with additive noise. The Annals of Probability, 38(2):532— 569, 2010. 

7. Jari Kaipio and Erkki Somersalo. Statistical and computational inverse problems, vol- 
ume 160. Springer, 2004. 

8. S. Lasanen. Discretizations of generalized random variables with applications to inverse 
problems. Ann. Acad. Sci. Fenn. Math. Diss., University of Oulu, 130, 2002. 

9. S. Lasanen. Measurements and infinite-dimensional statistical inverse theory. PAMM, 
7:1080101-1080102, 2007. 

10. Sari Lasanen. Non-gaussian statistical inverse problems, part i: Posterior distributions. 
Inverse Problems and Imaging, 6(2):215-266, 2012. 

11. Sari Lasanen. Non-gaussian statistical inverse problems, part ii: Posterior convergence 
for approximated unknowns. Inverse Problems and Imaging, 6(2):267-287, 2012. 

12. K.J.H. Law. Proposals which speed-up function spa<;e MCMC. arXiv preprint 
arXiv: 121 2.4767, 2012. 

13. Andrew C Lorenc. The potential of the ensemble kalman filter for nwpa comparison 
with 4d-vaj:. Quarterly Journal of the Royal Meteorological Society, 129(595) :3183-3203, 
2003. 

14. J. C. Mattingly. Ergodicity of 2D Navier-Stokes equations with random forcing and 
large viscosity. Commun. Math. Phys., 206:273-288, 1999. 

15. Andrew M Stuart. Inverse problems: a bayesian perspective. Acta Numerica, 19(1):451- 
559, 2010. 

16. L. Tierney. A note on Metropolis-Hastings kernels for general state spax;es. Ann. Appl. 

Probab., 8(1): 1-9, 1998. 

17. Sebastian J VoUmcr. Dimension-indepcuident MCMC sampling for elliptic inverse prob- 
lems with non-Gaussian priors. arXiv preprint arXiv:1302.2213, 2013. 

18. Dusanka Zupanski. A general weak constraint applicable to operational 4dvar data 
assimilation systems. Monthly Weather Review, 125(9):2274-2292, 1997. 



